o 



> 

00 



Adaptive and Recursive Time Relaxed Monte Carlo methods 

for rarefied gas dynamics 

Stefano Trazzi * Lorenzo Pareschi J Bernt Wennberg * 

October 3, 2008 



Abstract 

Recently a new class of Monte Carlo methods, called Time Relaxed Monte Carlo 
\ (TRMC), designed for the simulation of the Boltzmann equation close to fluid regimes 

y^, • have been introduced (TBJ- A generalized Wild sum expansion of the solution is at the 

basis of the simulation schemes. After a splitting of the equation the time discretization 
■ of the collision step is obtained from the Wild sum expansion of the solution by replacing 

high order terms in the expansion with the equilibrium Maxwellian distribution; in this 
way speed up of the methods close to fluid regimes is obtained by efficiently thermalizing 
particles close to the equilibrium state. In this work we present an improvement of such 
methods which allows to obtain an effective uniform accuracy in time without any 
restriction on the time step and subsequent increase of the computational cost. The 
main ingredient of the new algorithms is recursivity [18] . Several techniques can be 
used to truncate the recursive trees generated by the schemes without deteriorating 
£ — . the accuracy of the numerical solution. Techniques based on adaptive strategies are 

presented. Numerical results emphasize the gain of efficiency of the present simulation 
. schemes with respect to standard DSMC methods. 

o , 
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c5 1 1 Introduction 



Computations for rarefied gas dynamics (RGD) in engineering applications are most fre- 
quently performed using Monte Carlo methods. The Bird's method has been particularly 
successful for a wide range of applications [3j. One of the major drawbacks of these meth- 
ods is the difficulty to compute the simulation of rarefied gases that are close to the fluid 
dynamic limit, since in such regime the collisional time becomes very small. A nondimen- 
sional measure of the significance of collisions is given by the Knudsen number e, which is 



'University of Ferrara, Department of Mathematics, Via Machiavelli 35, 44f00 Ferrara Italy. E-mail: 
trazzi@dm.unif e . it 

' University of Ferrara, Department of Mathematics, Via Machiavelli 35, 44100 Ferrara Italy. E-mail: 
lorenzo .pareschi@unif e . it 

* Chalmers University of Technology, Department of Mathematics, SE-412 96 Goteborg, Sweden. E-mail: 
wennberg@chalmers . se 



1 



small in the fluid dynamic limit and large in the rarefied state. For small Knudsen numbers 
most Monte Carlo methods lose their efficiency because they are forced to operate on a very 
short time scale. The aim of this paper is to introduce a new Monte Carlo method that is 
robust in the fluid dynamic limit, by which we mean that it is accurate and efficient for a 
full range of Knudsen numbers. Alternative approaches to the derivation of efficient Monte 
Carlo methods for the simulation of the Boltzmann equation have been presented by several 
authors (see for example HUJ [HI 021 151 ES] and the references therein). Most of these 
approaches focus on the variance reduction of the method using information coming from 
the macroscopic scale (equilibrium states) at different levels. 

The starting point of the construction of this new family of schemes is the time relaxed 
discretization of the Boltzmann equation by the Wild sum [26] . Given a set of particles one 
tries to split it into subsets of particles according to the probabilities given by the coefficients 
in the sum. The remaining particles are sampled from a Maxwellian correspondent to the 
local equilibrium. This can be done with a recursive algorithm in an efficient way. As a 
result the method does not contain any time discretization error (similarly to Bird's method) 
on the contrary to Nanbu-Babovsky method or classical TRMC methods |144 116j. 

The goal of recursive TRMC (TRMC-R) methods is to construct simple and efficient 
numerical methods for the solution of the Boltzmann equation in regions with a large 
variation in the mean free path. As a consequence the TRMC-R methods have the following 
features: 

• for large Knudsen numbers, the TRMC-R methods behave as a Bird's method; 

• for intermediate Knudsen numbers the methods adapt the length of the collision trees 
in order to speed up the computation time without degradation of accuracy; 

• in the limit of the very small Knudsen number, the collision step replaces the distribu- 
tion function by a local Maxwellian with the same moments. The methods will behave 
as a stochastic kinetic scheme for the underlying Euler equations of gas dynamics |19| ; 




• mass, momentum, and energy are preserved. 

The paper is divided into 5 sections. After this introduction we will recall some basic 
facts about the Boltzmann equation and its fluid-dynamic limit. Next in section 3 we will 
discuss the problem of the time discretization of the Boltzmann equation and the basic 
notations of TRMC method. Section 4 is devoted to a detailed description of TRMC-R 
for Maxwell and Hard Sphere models. We present also adaptive techniques to improve the 
efficiency without degradation in accuracy. Finally in the last section we present detailed 
numerical tests both for a space homogeneous case and for a stationary shock wave, us- 
ing the Hard Sphere Model. The results show a marked improvement in the efficiency of 
computations given by TRMC-R over Bird's method. 

2 The Boltzmann Equation 

The Boltzmann equation describes the evolution of a continuum of particles by mean of 
three variables: the time t, the position x and velocity v of particles. 
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In this model, the density / = f(x,v,t) of particles follows the equation 

Of 1 , o 

-^ + v-V x f = -Q(f,f), x en cm. 3 , v em 3 , (l) 

supplemented with the initial condition 

f(x,v,t = 0) = f (x,v). (2) 

The function / can depend on other independent variables like an internal energy [7]. 

In ([T]) the parameter e > is called Knudsen number and it is proportional to the mean 
free path between collisions. The bilinear collisional operator Q(f, /) which describes the 
binary collisions between particles in a mono-atomic gas is given by 

Q(f,f)(v)= [ [ B(f(v')f(vi)-f(v)f(v*))dv*da, (3) 
Jr s Js 2 

where for simplicity the dependence of / on x and t has been omitted. 

In the previous expression a is a vector of the unitary sphere S 2 C R 3 . The collisional 

velocities (V,t>*) are associated to the velocities (v,v*) and to the parameter a by the 
relations 

v ' = ^ v + v * + M 0- )' v * = + v * + I 9 '* 7 )' ( 4 ) 

where q = v — is the relative velocity. 

The kernel B is a non negative function which characterizes the details of the binary in- 
teraction between particles. The classical Variable Hard Spheres model used for hypersonic 
flows in the upper-atmosphere is 

B(\q\,\q-a\)=K\q\ a ,0<a<l, 

where K is a positive constant. The case a = corresponds to a Maxwellian gas, while 
a = 1 is called a Hard Sphere Gas. 

Since / is a mass density in the phase space to obtain the density p = p(x, t) we have 
to integrate / in v 

P= ! fdv. (5) 
Similarly the gas velocity u is determined by 

pu = fvdv, (6) 
Jm 3 

and the gas temperature by 



3p 

Finally we define 



(v - uffdv. (7) 



E=\f \v\ 2 fdv=*-pT + \pu 2 (8) 

1 Jm3 2 Z 



the energy per volume density. 
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The collisional operator is such that the H-Theorem holds 



/ 



Q(f,f)log(f)dv<0. (9) 



This condition implies that each function / in equilibrium (i.e. Q(f, f) = 0) has locally the 
form of a Maxwellian distribution 

M( P ,u,T)(v) = ^f^exp [- L - W L ) , (10) 

where p, u, T are the density, the mean velocity and the gas temperature. 

Now, if we consider the Boltzmann equation ([I]) and multiply it for the elementary 
collisional invariants l,v, \v\ 2 and integrate in v we obtain 

Q(f,f)</>{v)dv = 0, cf>(v) = l,v,\v\ 2 , 

which correspond to conservation of mass, momentum and energy. If it is possible to invert 
the order of derivation - integration, we get 

| J <P(v)dv + £ JL J v^fdv = 0. (11) 

Unfortunately replacing <f> by the functions l,v, \v\ 2 we obtain a differential equation 
system which is not close since it involves higher order moments of the function f(x,v,t). 
Formally as e — > the function / is locally replaced by a Maxwellian. In this case it is 
possible to compute / from its moments using (I10p thus obtaining to leading order the 
closed Euler system of compressible gas dynamics 

i=i 

4(/™i) + Yl ~§7.( pUiU ^ + ~ir p = °' j = l > 2 ' 3 ( 13 ) 

i=i 1 i 
dE ^ d 

i=l 

where p = pT. 



3 Time Relaxed schemes 
3.1 Time discretizations 

The starting point is the usual first order splitting in time of (H|), which consists of solving 
separately a purely convective step (i.e., Q = in ([!])) and a collision step characterized by a 
space homogeneous Boltzmann equation (i.e., V x f = in §T§). Clearly, after this splitting, 
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almost all the main difficulties are contained in the collision step. For this reason, in what 
follows we will fix our attention on the time discretization of the homogeneous Boltzmann 
equation 

!=>,/)■ m 

As proposed in |8j, a general idea for deriving robust numerical schemes, that is, schemes 
that are unconditionally stable and preserve the asymptotic of the fluid-dynamic limit, for 
a nonlinear equation like (|15p . is to replace high order terms of a suitable well-posed power 
series expansion by the local equilibrium. Here we will briefly recall the schemes presented 
in [8]. 

3.2 Derivation 

Let us consider a differential system of the type 

% = - £ [P(fJ)-»f}, (16) 

with the same initial condition ([2]), and where /i ^ is a constant and P a bilinear operator. 
Let us replace the time variable t and the function / = f{v,t) using the equations 

r = (1 - e-^/ £ ), F(v, r) = f(v, t)e^ . (17) 

Then F is easily shown to satisfy 

dF 1 

^ = -P(F,F) (18) 

with F(v,t = 0) = f (v). 

Now, the solution to the Cauchy problem for (|18p can be sought in the form of a power 
series 

oo 

F(v,t) = Y^r k f k (v), f k=0 (v) = f (v), (19) 

A:=0 

where the functions are given by the recurrence formula 

fk+l(v) = -r^- 1 Y / ~P(fh,fk-h), k = 0,1, (20) 

Making use of the original variables, we obtain the following formal representation of 
the solution to the Cauchy problem (fT5|) : 

oo 

f(v t t) = e-^ £ ]T(1 - e-^) k fk(v). (21) 

fc=0 

The method was originally developed by Wild |26} [6] to solve the Boltzmann equation 
for Maxwellian molecules. 

From this representation, a class of numerical schemes can be naturally derived. 
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In [S], the following class of numerical schemes, based on a suitable truncation for m > 1 
of dH|), has been constructed: 

m 

f n+1 (v) = e -» At/£ - e^ A * /e ) fc /r(^) + (1 - e-^ A * /e ) m+1 M(t;), (22) 

fc=0 

where / n = f(nAt) and At is a small time interval. The quantity M (referred to as the 
local Maxwellian associated with /) is the asymptotic stationary solution of the equation. 

It can be shown that the schemes obtained in this way are of order m in time. Further- 
more, these schemes satisfy the following properties [8]. 

(i) Conservation. 

If P(f, g) is a nonnegative bilinear operator such that there exist some functions <f>{v) 
with the following property, 

/ P(f,f)4>(v)dv = l i [ f<p{v)dv, (23) 

and the initial condition f° is a nonnegative function, then f n+1 is nonnegative for 
any fj,At/e and satisfies 

/ f n+1 4>(v)dv= I f n Hv)dv. (24) 

JRZ JR? 

(ii) Asymptotic preservation (AP). 
For any m > 1, we have 

lim f n+1 = M(v). (25) 

/j,At/e— >oo 

In the case of the Boltzmann equation, with a collision kernel bounded by a taking 
P{fi f) = Qifi f) + A 4 /] with fj, > Airpa the schemes guarantee the conservation of mass, 
momentum and energy (by the first property) and the correct solution near the fluid limit 
(i.e. e -»• 0). 

4 TRMC-R methods 

4.1 Maxwellian case 

In order to simplify the derivation of the Recursive Time Relaxed Schemes, we recall some 
basic facts on the algorithm for the simple case of constant cross sections (Maxwellian 
molecules) as proposed in [16] 

First we note that in the case of Maxwell molecules the Wild sum has a clear prob- 
abilistic interpretation. If /(At) is the velocity distribution of particle at time At, then 
taking a particle at random from this distribution it might happen that this particle has 
not collided one single time. The distribution given this is just /o and the probability 
of this event is exp(— /xAt/e). In the same way f\ is the velocity distribution for par- 
ticles which have been involved in exactly one collision, and the probability of that is 
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(1 — exp(— /iAt/e))(exp(— fxAt/e)). At least f m is the velocity distribution given that ex- 
actly m + 1 particles have been involved in their collision history back to the initial time. 
To be able to find a sample of f m , we must assume that the densities /fc,0 < m — 1 are 
all already known. Off course the only one of these that is really known is /o, the initial 
distribution. 

A sample of f m can be determined in a recursive way. In order to understand how the 
algorithm works, let's consider for example the sampling from a starting density functions 
fk, k = 1,2,3. From (|20p fx is obtained from the collisional operator P(fo,fo), f2 from a 
combination of P(fo, fx) and P(/i,/o) with the same probability weight, where the terms 
/i are constructed as seen before. Similarly /2 is a combination of P(fo, P(fo, /o)) and 
P(P(/o,/o),/o). 

In the same way it is possible to create the higher terms of the Wild's Sum. 
A simple recursive Monte Carlo algorithm is the following 



Algorithm 4.1 (Recursive sampling for Maxwell molecules) 

1. choose n from a geometric distribution with parameter 
t = 1 — exp(— /iAt/e) 

2. take a sample from the distribution with density f n 

o if n = 0, take v from the initial density /o 
o else proceed as follows 

- choose k G {0, 1, . . . , n — 1} with equal probability 

- take a sample Vi from the density 

and Vj from the density f n -k-i as in step 2 

- perform the collision between Vi and Vj , obtaining v[ and v'j 

- then d- and v'j are distributed according to the density f n 



The post collisional velocities are computed through relations 

, Vi+Vj \vi — Vj\ , Vi + Vj \Vi-Vj\ 

Vi = — + i a > V J = 2^ ^ ' 

where u) is chosen uniformly in the unit sphere, according to 



cos <f> sin 

sin0sin6l ), 9 = arccos(2^ - 1), = 2vr£ 2 , (27) 
cos 9 

and £1,^2 are uniformly distributed random variables in [0, 1]. 

It is useful to use a representation of the collision process through the collision trees, 
sometimes called Mc Kean graphs (Figure [TJ. 

It is now clear that the velocity distribution fk can be drawn from the starting distri- 
bution /o, choosing the different collisional trees with the same probability in mean. Off 
course two particles are produced in every collisional event, and only one of these is used 
to complete the collisional process. So it is natural to store the particle not used in a list, 
for its direct utilization in a future collisional process. In this way we can increase the 
efficiency of the method (if a particle sampled from a velocity density fk already exist it 
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Figure 1: Mc Kean graphs for f±, fi and fj, 



will not be necessary to obtain it by the complete collisional process) and guarantee the 
exact conservation of moments. 

In order to do this we first split a set of N particles into m collision sets Ni, i = 0, . . . , m 
where each set iVj characterizes the number of particles that will undergo i collisions in a 
time step At. All particles having more then m collisions in the time step are thermalized 
(i.e. replaced with an equilibrium particle taken from the local Maxwellian), this number 
is denoted by N m+ \. Accordingly to this we have the following algorithm where N and m 
are given. 



Algorithm 4.2 (Splitting particles into collision sets) 

1. set t = 1 — exp(— /zAt/e) 

n = 0, A = I, u Q = (1 - r), N Q = N, N c = 

2. repeat 

- n = n + 1, 

- N n = 7V„_i - N n -! 

- L> n = (1 - t)t" 

- A n = A„_i — UJn-l 

- N n = Iround(*^) 
-N C = N C + N n 

until N n > and n < m 

3. if (n — m) 

- N rn+1 = N — N c — N 
else 

- m = n, N m+ i = 
end if 
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Here, by Iround(x), we denote a suitable integer rounding of a positive real number x. 
In our algorithm, we choose 

, _ f [x] with probability [x] + 1 — x, 
\ [x] + 1 with probability x — [x], 

where [x] denotes the integer part of x. 

Note that, since we have a finite number of particles, a consequence of the above splitting 
into collision sets is that the maximum possible length of a collision process is fixed by the 
initial number of particles N. The recursive collision algorithm for Maxwell molecules 
where thermalization occurs accordingly to a TR discretization of order m is given here. To 
achieve exact conservation of momentum and energy and a better computational efficiency, 
the algorithm uses counters to keep track of the number of particles stored in memory 
with a collision history of i collisions in a time step At and not yet used in the simulation. 



Algorithm 4.3 (TRMC-R for Maxwell molecules) 

1. compute the initial velocity of the particles, {v® ,i = 1, . . . , iV} 
by sampling them from the initial velocity fo 

2. split particles into collision sets as in algorithm 

3. set counters c n — for n — 1, . . . , m + 1 
4- for n = m, . . . , 1 

take N n samples from the distribution with density f n , according to 
o repeat 

- choose k 6 {0, 1, . . . , n — 1} with equal probability 

- if k = take Vi from the initial density fo 

- else choose Vi from the density fk 

if Ck > use a stored particle with a random choice 

set Ck = Ck — 1 and Nk — Nk + 1 

else sample Vi and v* from fk (recursively) 

v* is stored and then set Ck = c/. + 1. Nk = Nk — 1 

- if n — k — 1 = take Vj from the initial density fo 

- else choose Vi from the density f n -k-i 

i/cn-fc-i > use a stored particle with a random choice 
set c n -k-\ = c n -k-i - 1 and N n -k-i = N n -k-i + 1 
else sample Vi and v* from f n -k-\ (recursively) 
v* is stored and then set c n -k-i = c n -k-i + 1, 
N n -k-i = N n -k-i — 1 

- perform the collision between Vi and Vj as in DSMC 

- v'i and v'a are random variables distributed 
according to the density f n 

- set N n = N n - 2 
until (N n > 0) 

end for 

5. sample N m+ i particles from the local Maxwellian 



9 



4.2 VHS collision kernels 

The algorithm described for Maxwellian molecules can be extended to more general collision 
kernels by using dummy collisions and acceptance-rejection technique. This approach is 
equivalent to sample the post collisional velocity according to P(f, /) / \i, where \x = Ana 
and a is an upper bound of the scattering cross section for the given set of particles. 

The upper bound a should be chosen as small as possible, to avoid inefficient rejection, 
and it should be computed fast. 

An optimal bound can be derived taking a as 

a = maxa(\vi — Vj\). (28) 

Vi,Vj 

However this computation would be too expensive since it would require an 0(N 2 ) opera- 
tions. An upper bound of a can be obtained by taking a = a(2Av), 

Av = max \vi — v\, v = V^Uj/JV. 

i 

In the VHS case the algorithm should be modified as follows 

Algorithm 4.4 (TRMC-R for VHS molecules) 

1. compute the initial velocity of the particles, {v® ,i — 1, . . . , iV} 
by sampling them from the initial velocity f 

2. split particles into collision sets as in algorithm \4-<fy 

3. set counters c n = for n = 1, . . . , m + 1 

4- compute an upper bound a of o~ij = o~(vi,Vj) 

5. for n = m, . . . , 1 

take N n samples from the distribution with density f n , according to 
point 4 of algorithm \4-3\ where the dummy collision is performed as 
in DSMC if a£± < Oij, with £i uniformly distributed random variable in [0, 1] 

6. sample N m +i particles from the local Maxwellian 

The algorithm is exactly conservative if combined with a suitable scheme for sampling 
a set of particles with prescribed momentum and energy from the Maxwellian, as proposed 
in [19] or by the authors in [17] . 

Remark 4.1 

• In the case of m = 1,2 the TRMC-R method corresponds to the first and second 
order TRMC methods presented in fToJ for the simplest possible choice of the weights. 
Off course here we are aiming at using much larger values of m 3> 1 for which a 
direct extension of the algorithms presented in ]16^ is not feasible. Larger values of m 
produce higher accurate results and then allow a larger time step when compared to 

m\ns- 

• Since the collision process changes the distribution function, it is important to choose 
a correct upper bound a of aij in order to avoid discarding all the collisional trees 
computed in case o~ij > a. As an alternative one can update the upper bound itself, 
after each collision as in ]lb^ \17f . 
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4.3 Adaptive technique for TRMC-R scheme 

In practical simulations the number m can be very large, depending on the Knudsen num- 
ber and on the number of test particles. Clearly small values of m make the algorithm 
faster, because the collision process is replaced by the projection to the local Maxwellian 
equilibrium, but far from the fluid regime keeping m too small can produce less accurate 
results. In practice, a maximum allowed value m max of m is fixed at the beginning of the 
calculations; m max represents the maximum depth of a collision tree. The main problem 
is to choose the right m max , in order to have the best combination between efficiency and 
accuracy. The idea we develop is to use an adaptive technique to choose the right maxi- 
mum depths of the collision trees, based on evaluating the distance of the solution from the 
equilibrium through a suitable indicator. This can be performed measuring the variation 
of some macroscopic variables such as the fourth order moment or the components of the 
shear stress tensor. 

Let S be the macroscopic variable selected according to the particular physical problem. 
Then define the quantity 

|gn+l,m m ax _ gnj 



\s r - 



that represents the relative variation at time step n + 1 of the macroscopic variable S 
computed with the solution obtained using m max as maximum depth of the collision trees. 
If we fix an interval [Si, 8%], Q < 5\ < 82 then we can apply the following criteria in order to 
accept or discard the solution at time step n + 1 

• if Ei < 5i the solution is accepted and m max = m max /2 in the next time step n + 2; 

• if £1 < Ei < 82 the solution is accepted and rra max is unchanged in the next time step 
n + 2; 

• if Ei > 82 the solution is discarded and the calculation is performed again using 

Timax = 2f7l max . 

Off course it is possible to use other techniques in order to evaluate the distance from 
equilibrium |24j . 

The algorithm works with optimal efficiency if the collisions computed with the "wrong" 
m max are kept and reused with 2m max . 

This can be done by starting with m = m m i n and observing that if 



fn+l,m = T - T ) k f k + (1 - r) m+1 M 



k=0 

then 

2m 

jn+l,2m = Jn+l,m _j_ j- ^ (j _ T )fcj fe + [(1 _ r )2m+l _ ^ _ T )m+lj M 
fc=m+l 

So if the test for the adaptive strategy impose to discard the solution we proceed as 
follows 
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Definition 


left tree 


right tree 


([29]) 


7 


7 


m 


3 


2 


m 


3 


43/16 



Table 1: Length of the collision trees 



• the collisions computed with m are kept; 

• the fraction (1 — r) m+1 M is discarded; 

• the fraction ^f= m+ i(l ~~ T ) k fk is computed by the recursive collision process; 

• the fraction [(1 — r) 2m+1 — (1 — r) m+1 ] is sampled by a Maxwellian. 

We observe that for estimation purposes the sampling from Maxwellian can be substitute 
by the analytical computation of local Maxwellian. This contributes to a better efficiency 
of the adaptive method. 

4.4 Truncation of the collision trees for TRMC-R scheme 

Off course several definitions of the length of trees are possible; the simplest one, that does 
not care about the shape of the collision trees that generate a particle from the density 
function is the one provided by 

L(k = h + j + 1) = k. (29) 

The length corresponds to the coefficient of the density function by which we sample the 
particle. The idea is to sample directly from the local Maxwellian if, for the collision 
process, L(k) > m max , m max fixed, otherwise the whole tree is kept and the collisions 
are performed. This simple definition has been used into the algorithms described before. 
Different definition of length L can be done using recursivity as 

L(k = h + j + l) = l + min{L(/i),L(j)}, (30) 
L{k = h + j + l) = 1 + mean{L(/i),L(j')}. (31) 

These two last definitions can be related to the concept of "well balanced and not well 
balanced trees" (see Figure [2]), i.e. if L(k) > m max at the end of the collision process we can 
imagine that the particles would be more thermalized with respect to the ones that come 
from a tree where L(k) < m max . This off course is not true for definition (|29p . 

In table ([T|) we show the values of the length of the collision trees reported in Figure ([2]) 
using the three different definitions. 

The implementation of such strategy inside the recursive algorithm needs a modification 
of the first formulation given in section 14.21 because it is necessary to evaluate the length 
of collision trees without performing collisions. The idea is to write into a list the collision 
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Figure 2: Well balanced (left) and not well balanced tree (right) 

process (performed in a recursive way), to evaluate the length and to perform the collision 
using the collision process stored into the list. 

The following algorithm shows how we can store the collision tree (correspondent to a 
process to sample from the density function /*.) into a list named tree and how we can 
evaluate his length by using the recursive definition seen above. 

Algorithm 4.5 (Storage and length of collision trees) 

1. assign to a variable indx the value indx = 

2. put path(indx) = k 

3. if k = return 

else proceed as follows 

- choose j £ {0, 1, . . . , k — 1} with equal probability 

- put h = k — j — 1 

- increment the index: indx = indx + 1 

- store j in the list: path(indx) = j 

- increment the index: indx = indx + 1 

- store h in the list: path(indx) = h 

- return (1+ 

+ mh\{{repeat from step 3 with k — j), (repeat from step 3 with k = h)} 



In this way we obtain a list path which contains the structure of the collision tree and its 
relative length using definition (|3ip . Off course the same algorithm applies also to definition 

dSU). 

Now we can finally modify the algorithm 14. 4l at the point 5 in order to use the previously 
created collision trees. 

Algorithm 4.6 (TRMC-WB for VHS molecules) 

1. compute the initial velocity of the particles, i = 1, . . . , iV} 
by sampling them from the initial velocity fo 

2. split particles into collision sets as in algorithm \4-.2fy 

3. set counters c n — for n = 1, . . . , m + 1 
4 ■ compute an upper bound a of o~ij 

5. for n = m, . . . , 1 
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take N n samples from the distribution with density f n , according to 
repeat 

- compute the list path and its length L as in algorithm \4-5\ 

- ifL> ^max 

set N rn+ i = N rn+ i + 1 and N n = N n — 1 

else 

set indx = 0, k = path(indx) 
5.1. if k = sample from the initial density f Q 
else indx — indx + 1 

set j = path(indx) 

set indx = indx + 1 

set h = path{indx) 

if Cj > use a stored particle with a random choice 

set Cj = Cj — 1, Nj = Nj + 1 and indx — indx + 2j 

else sample Wj a-nd v* from fj ( recursively from 5. 1 
by setting k = j) 

v* is stored and then set Cj = Cj + 1, Nj — Nj — I 

end if 

if Ch > use a stored particle with a random choice 

set Ch = Ch — 1, Nh = Nfr + 1 and indx = indx + 2h 

else sample Vj and v* from fh (recursively from 5.1 
by setting k = h) 

v* is stored and then set Cj = Cj + 1, Nj = Nj — 1 

end if 

Perform the dummy collision between Vi and Vj as in DSMC. 
end if 

set N n =N n -2 

end if 
until (N n > 0) 
end for 

6. sample N m+ i particles from the local Maxwellian 



Remark 4.2 In the Variable Hard Sphere case the calculation of algorithm \4-5\ does not 
represent the effective length of the real collision trees, due to the dummy collisions which 
can occur during the total collision process. Thus the collisional trees performed in VHS 
model are usually longer than the ones for Maxwell molecules. 

5 Numerical tests 

We present some numerical tests, both in space homogeneous and space non homogeneous 
situations using the Hard Sphere model. The solution obtained by the Recursive Time 
Relaxed scheme combined with the adaptive strategy, is compared with the one obtained 
by the classic Bird's algorithm. 

To simplify notation we will use TRMC-R for the basic Recursive Time Relaxed Monte 
Carlo Scheme defined by Algorithm 4.4, TRMC-RAD for same scheme improved by the 
adaptivity strategy of section 4.3 (based on the shear stress tensor as equilibrium indicator) 
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and TRMC-WB for the Recursive Time Relaxed Monte Carlo Scheme defined by Algorithm 
4.6 (based on the well balanced truncation of the collision trees). With BIRD we refer to 
the classic DSMC scheme [3]. 

All macroscopic quantities have been considered in non dimensional form. As efficiency 
indicator we have considered the total number of collisions performed in the simulation. 
Sampling two particles from the Maxwellian has been counted as one collision. Note that 
due to the structure of the recursive algorithm in principle we may have additional compu- 
tational requirements in terms of storage of particles during the evaluation of collision trees 
and storage and calculation of the collision trees if we want to implement a shape-depending 
truncation of the trees. However particle storage along the collision trees does not require 
any additional memory. In fact to achieve conservations each particle is used only one time 
in a collision tree and so it is the particle itself which is stored in the list. On the other 
hand the additional use of memory for the storage of the computational trees involves only 
one single collision history and thus the impact in the overall calculation is negligible. 

5.1 Space homogeneous tests 

We consider the sum of two Maxwellian as initial data. The solutions has been obtained in 
one single run using 5 x 10 4 particles and choosing a collision time step At/e equal to one. 
The reference solution, called in the plot as "exact" has been obtained by a large number 
of averages on Bird's scheme. 

We compare the results for the fourth order moment M4 and the component P xx of the 
shear stress tensor 

M 4 = f fvUv, P xx = [ f( Vl - Ul ) 2 dv. 

JR3 J R 3 

Figure U] shows the case of TRMC-R without any limit for the maximum depth of 
collision trees. The agreement with respect to Bird's scheme is very good, and if we look at 
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Figure 5: Relative error for the evolution in time of the stress tensor component P xx (left) 
and number of collisions (right) for TRMC-R and Bird's method. 



the relative errors we can observe that they are comparable. In Figure [5] we plot the relative 
error for the shear stress tensor component P xx and the number of collisions necessary to 
perform the simulations. We can consider the number of collisions as an index of efficiency 
in terms of computational time. A single sample from a Maxwellian is counted as half of a 
collision. As expected we had no gain in efficiency because with this scheme we took into 
account all possible collision processes. 

In the next case (Figure [6] and [7]) we perform the same simulation using TRMC-RAD, 
starting with m max = 2 and taking 5i = 0.005 and 82 = 0.01. We preserve accuracy with 
respect to DSCM solution, but we have obtained a strong reduction of the computational 
time because a bigger fraction of particles is sampled directly from the Maxwellian without 
going through the whole collision tree. 

We observe in Figure [81 where we report the maximum depth of the collision trees 
during the calculation, that the simulation in the first time steps has been performed several 

times, Using ctt 6cich tim6 TTimax 

= 2m max . As already mentioned let us point out that these 
calculations do not affect the total computational cost because all "old collisions" are kept 
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Figure 6: Evolution in Time of 4 order moment (left) and relative error (right) for TRMC- 
RAD and Bird's method. 




Figure 7: Relative error for the evolution in time of the stress tensor component P xx (left) 
and number of collisions (right) for TRMC-RAD and Bird's method. 

and reused in order to complete the new collision process. 

The last homogeneous test deals with the shape depending strategy based on truncation 
of collision trees by using the definition L(k — /j.-t-j-t-l) — l-|-min{Zv(/i), Zv(j)} and ?7imax — 5. 
Similar results can be obtained using the definition L(k = h+j+1) = l+mean{L(/i), L(j)}, 
we omit them for brevity. We have, also in this case, a gain in efficiency, preserving the 
accuracy in time (see Figures [9] and fTUl) . 

The TRMC-WB scheme shows a constant gain of computational cost with respect to 
Bird's scheme, approximatively the 7% during the whole simulation, while the TRMC-RAD 
achieves the maximum gain (about 50%) at the end of the calculation, with an average gain 
of 20% in the total simulation time. The combination of TRMC-WB and TRMC-RAD is 
under investigation and will be presented elsewhere. 



17 




2 4 6 8 10 12 14 



Figure 8: Maximum depth m max of a collision tree in time for TRMC-RAD 
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Figure 9: Evolution in Time of 4 order moment (left) and relative error (right) for TRMC- 
WB and Bird's method. 

5.2 Stationary shock 

Next we consider a non homogeneous stationary shock problem. The boundary conditions 
have been assigned at the left and right side accordingly to the Rankine-Hugoniot relations 

Plul = Prur, 

PLU 2 L +PL = PRU R + PR, 

u L {E L + p L ) = u R (E R + p R ). 
The values used in the simulation are 

• M aL = 3 (Mach Number of incoming flux) 

• Tl = 1 (Temperature of incoming flux) 

• u XL = —My/^yT£, u yL = 0, u ZL =07 = 5/3 (Mean velocity of incoming flux) 

• pl = 1 (Total mass) 
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Figure 10: Relative error for the evolution in time of stress tensor component P xx (left) and 
number of collisions (right) for TRMC-WB and Bird's method. 

• e = 1, e = 0.1 and e = 0.001. 

The numerical solution has been obtained using TRMC-RAD and Bird's method with 50 
spatial cells and 1000 particles in each cell. Reference 'Exact' solution has been performed 
by Bird's standard DSMC method using 50 spatial cells and 3000 particles in each cell. A 
detailed analysis on the effect of the number of particles per cell in TRMC methods has 
been performed in [211 122] . 

In order to increase accuracy, for t large enough, averages of the solution have been com- 
puted. The results for the temperature shock are presented in Figures [TTlfTBl As expected 
there is a good agreement between TRMC-RAD and Bird's method and the relative errors 
are essentially comparable. 

Note that there is an evident gain in efficiency of TRMC-RAD against Bird's method 
without losing accuracy, especially near fluid regime. Looking at the rarefied case e = 1 we 
obtain the same computational cost, while we have a gain close to 10% for the intermediate 
test case e = 0.1 that increases up to the 86% close to the fluid regime for e = 0.001. 

6 Conclusion 

We have presented recursive Monte Carlo methods which are suitable for the numerical 
simulation of the Boltzmann equation for a wide range of Knudsen numbers. These recursive 
TRMC methods minimize the effects of time discretization error and over-relaxation due 
to the choice of the upper bound of the cross-section that were present in the previous 
versions of TRMC. Of paramount importance to increase the efficiency of the methods is 
the use of suitable truncation strategies of the collisional trees, such as adaptive truncation 
based on macroscopic quantities or well balanced truncation based on the trees properties. 
The resulting schemes are very promising and show that a considerable gain in efficiency 
can be obtained without degradation in accuracy. However additional test cases must be 
performed in order to validate the schemes. A combination of the present schemes with the 
hybrid strategy proposed in |15j is actually under investigation. 



19 




Comparison on number o! collisions tore = 1 
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Figure 12: Rarefied regime (e = 1). Number of collisions for TRMC-RAD and Bird's 
method (left) and maximum depth m max in each cell for TRMC-RAD (right). 




Figure 13: Intermediate regime (e = 0.1). Temperature (left) and Relative Error (right) for 
TRMC-RAD and Bird's method. 
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Figure 14: Intermediate regime (e = 0.1). Number of collisions for TRMC-RAD and Bird's 
method (left) and maximum depth m max in each cell for TRMC-RAD (right). 

Temperature fore = 0.001 Relative Error on Temperature fort = 0.001 




Figure 15: Fluid regime (e = 0.001). Temperature (left) and Relative Error (right) for 
TRMC-RAD and Bird's method. 



Comparison on number of collisions for e = 0.001 
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Figure 16: Fluid regime (e = 0.001). Number of collisions for TRMC-RAD and Bird's 
method (left) and maximum depth m max in each cell for TRMC-RAD (right). 
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